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FFTs and 
Beyond 



Digital signal 
processors are often 
used for doing 
spectral analysis of 
incoming waveforms, 
but the old standby 
FFT can have its 
downside in the real 
world. Find out what 
to watch for in your 
next design. 
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he analysis of a 
^^^B ^^^r signal based on its 

^^■^^ frequency content is 
^^commonly referred to as 
spectral analysis. Although the 
mathematical basis for this operation, 
the Fourier Transform, has been 
known for many years, it was the 
introduction of the Fast Fourier 
Transform (FFT) algorithm which 
made spectral analysis a practical 
reality. 

Implementing the FFT in personal 
computers and embedded DSP systems 
has offered an efficient and economical 
application of Fourier techniques to a 
wide variety of measurement and 
analysis tasks. Moreover, because the 
FFT has been found to be so valuable 
in applications such as medical signal 



processing, radar, and telecommunica 
tions, DSP chips are often designed to 
implement the FFT with the greatest 
efficiency. 

In most instances, the powerful 
Fourier techniques, used in modems, 
fax machines, and CT or ultrasound 
scanners, are hidden from the user, 
who doesn't have to worry about theii 
mathematical implications. In other 
cases, however, human interpreters 
must make diagnostic decisions based 
on frequency-domain representations 
of data processed through Fourier 
transforms. 

For example, many digital storage 
oscilloscopes offer the user the option 
of converting time-domain signals into 
the frequency domain through the use 
of the FFT which runs on an embedded 
DSP and displays results directly on 
screen. It is also common for scientists 
and engineers to write short tt 1 -based 
routines to display a spectral represen- 
tation of experimental data acquired 
by a personal computer. It is in these 
cases where the unwary may fall into 
one of the many traps that the FFTs 
conceal. 

FFT users often forget that real- 
world signals are seldom periodic, free 
of noise and distortion, and that signal 
and noise statistics play an important 
role in their analysis. Because of these 
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Figure 1— A purely sinusoidal signal (a) has a single impulse as ils true spectrum (b). However, the signal is viewed 
by FFT through a Unite window (c), and it is assumed that this record is repealed beyond the FFTs window (d). This 
leads to leakage of the main lobe to sidebbes in the spectral estimate (e). 
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"problem" factors, the FFTs and other 
methods can only provide estimates of 
the actual spectrum of signals. The 
results require competent interpreta- 
tion by the user for correct analysis. 

In this article, I will explain the 
common pitfalls in the use of the FFT 
and how to avoid them. After exposing 
some of the inherent problems which 
make the FFT unsuitable for high- 
resolution applications, I'll present 
more powerful spectral estimation 
methods, which cope with the funda- 
mental shortcomings of the FFT, and 
describe typical applications for these 
methods. 

FFTs AND THE POWER 
SPECTRAL DENSITY 

Using a typical data acquisition 
setup, a signal is sampled at a fixed 
rate of 

, f samples ! 
I second I 

which yields discrete data samples x , 
x,, x N These N samples are then 
equally spaced by the discrete sam- 
pling period Atfseconds] = The 
discrete Fourier transform (DFT) 
represents the time-domain data with 
N-spaced samples in the frequency 
domain X,„ X,, X N , through: 



X(f)»At Z x„e J "' , " 4t 



(M 



where the frequency /[Hz] is defined 
over the interval -Vm <f< Vz*. The 
FFT efficiently evaluates this expres- 
sion at a discrete set of N frequencies 
spaced equally by Af\Hz\ - '/«*. 

In its most simple form, the 
energy-spectral-density estimate of the 
time-domain data is given by the 
squared modulus of this data's FFT, 
and the power spectral density (PSD] 
estimate P(f) at every discrete fre- 
quency / is obtained by dividing the 
latter by the time interval NA1: 

HU-'^f ;m-0, 1 K-J (2) 

where /jHz) - mAf. In a case which 
uses real data (this is the norm when 
sampling from real-world signals), the 
PSD for negative frequencies is 
symmetrical to the PSD for positive 
frequencies, making only half of the 
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Table ^Typical window functions tor use with the FFT in spectral estimation. Windows are N points long and are 
assumed here to be symmetric around n = 0. 



PSD useful. However at times, it may 
be necessary to compute PSD for 
complex data where relevant results 
are obtained for both positive and 
negative frequencies. 

Although obtaining the PSD 
seems to be as simple as computing 
the FFT and obtaining the square 
modulus of the results, it must be 
noted that, because the data set 
employed to obtain the Fourier 
transform is a limited record of the 
actual data series, the PSD obtained is 
only an estimate of the true PSD. 
Moreover, as will be seen later, 
meaningless spectral estimates may be 
obtained by using Equation |2) without 
performing some kind of statistical 
averaging of the PSD. 

PITFALLS OF THE FFT 

When sampling a continuous 
signal, information may he lost 
because no data is available between 
the sample points. As the sampling 
rate is increased, a larger portion of the 
information is made available. Accord- 
ing to Nyquist's theorem, to correctly 
sample a waveform, the sampling rate 
must be at least twice that of the 



highest-frequency component of the 
waveform. Disregarding this rule will 
result in aliasing — a process in which 
signal components of frequency higher 
than half the sampling rate appear as 
components with a frequency equal to 
the difference between the actual 
frequency of the component and the 
sampling rate. 

Because aliased components 
cannot be distinguished from real 
signals after sampling, aliasing is not 
just a minor source of error. It is 
therefore of extreme importance that 
antialiasing filters with very high roll- 
off be used for all serious spectral 
analysis. 

Beyond appropriate sampling 
practices, the FFT still exposes other 
inherent traps which can potentially 
prevent analysis of a signal. The most 
important problems include leakage 
and the picket-fence effect. 

Leakage is caused by the fact that 
the FFT works on a short portion of 
the signal, a phenomenon called 
windowing, because the FFT can only 
see the portion of the signal that falls 
within its sampling "window," after 
which it assumes that windowed data 



The Computer Applications Journal Issue #52 November 1 994 



21 



repeats itself indefinitely. However, as 
shown in Figure 1, this assumption is 
only seldom eoneel. In most eases, the 
FFT analyzes a distorted version of the 
signal that emu. mis discontinuities 
resulting from appending windowed 
data to their duplicates. In PSD, these 
discontinuities appear as a leakage of 
the energy's real frequency compo- 
nents into sidelobes which show up on 
either side of a peak. 

The second problem, called the 
picket-fence effect or scalloping, is 
inherently related to the discrete 
nature of the DFT. That is, the DFT 
calculates the frequency content of a 
signal at very well-defined discrete 
points in the frequency domain rather 
than producing a continuous spec- 
trum. In a perfect system, if a certain 
component of the signal had a fre- 
quency falling between the discrete 
frequencies computed by the DFT, this 
component would not appear in the 
estimated PSD. 

To visualize this problem, suppose 
that an ideal signal is sampled at a rate 
of 2048 Hz and processed through a 




-P(0 



Figure 2— In one implementation ol a parametric high-resolution spectral estimator, the coefficients a t , a, a p of 
anAR tiller are determined from input data through Marple 's algorithm. The transfer lunction of the filter is then 
evaluated efficiently by FFT, resulting in a high-resolution estimate of the input data's PSD. 



256-point FFT. There would be a 
spectral channel every 4 Hz [at DC, 4 
Hz, 8 Hz, 12 Hz, etc.). Suppose now 
that the signal being analyzed is a pure 
sinusoidal with a frequency of 10 Hz. 
In a perfect system, this signal would 
not appear in the PSD because it falls 
between two discrete frequency 
channels — much like the case of a 



p 



picket fence which allows us to see 
detail in the scene behind it only if the 
details happen to fall within a slot 
between the boards. 

In reality, however, because the 
FFT produces slightly overlapping 
"bins" of finite bandwidth, compo- 
nents with frequencies that fall 
between the theoretical discrete lines 
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arc distributed among adjacent bins, 
but at reduced magnitudes. This 
attenuation is the actual pickct-icncc 
or scalloping error. Both of these 
problems arc somewhat corrected by 
the use of an appropriate window. 

So far, all samples presented to the 
FFT have been considered equal, 
which means that a weight of one has 
been implicitly applied to all samples. 
The samples outside of the FFT's scope 
are not considered, and thus their 
effective weight is zero, resulting in a 
rectangular-shaped window. This 
ultimately leads to the discontinuities 
that cause leakage. 

A number of windows have been 
devised which reduce the amplitude of 
the samples at the edges of the 
window while increasing the relevance 
of samples towards its center. By doing 
so, these windows reduce the disconti- 
nuity to zero, thus lowering the 
amplitude of the sidelobes that 
surround a peak in the PSD. In 
addition, the use of a nonrectangular 
window increases the bandwidth of 
each bin, which results in a decreased 
scalloping error. 

Some typical window functions 
and their characteristics are presented 
in Table 1. In essence, these functions 

produce N weights w„, w, w N , 

which arc "weighted" (multiplied) 
one-to-one with their corresponding 
data samples x ff x,, x N , before 
subjecting them to the FFT: 

Xffj.AtZ' v^x„e-^' M (3) 

n.O 

Reduced resolution is the price 
paid for a reduction in leakage and 
scalloping through the use of a 
nonrectangular window. In fact, if it is 
necessary to view two closely spaced 
peaks, the rectangular window's 
narrow main lobe lets the user obtain 
analysis results, which report the 
existence of these closely spaced 
components. Any of the other win- 
dows would end up fusing these two 
peaks into a single smooth crest. 

The use of a rectangular window 
is also appropriate for the analysis of 
transients. In these cases, a zero signal 
usually precedes and succeeds the 
transient. Thus, if the FFT is forced to 
look at the complete data record for 
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Figure 3a— Positive frequency theoretical spectrum of Mamie's 64-point complex data test set. This spectrum 
includes features that are well suited for evaluating spectral estimators. 



the transient, no artificial discontin- As you see, there is no single 

uities are introduced, and full resolu- window which outperforms all others 
tion can be obtained without leakage. in every respect, and it is safe to say 



PSD estimates 
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Figure 3b— Unlike the very well-established features of the theoretical spectrum shown in Figure 3a, spectral 
estimates distort the PSD according to their inherent assumptions. The spectrum is estimated using three different 
methods: i) zero-padded FFT pehodogram (green line), ii) Welch's method (black line. N = 32,S= 16), and Hi) 
Mamie's method (red line, p = 15). 
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that selecting the appropriate window 
for a specific application is more of an 
art than an exact science. 

Another solution comes in handy 
when the signal rides on a relatively 
high DC level or on a strong sinusoidal 
signal. In these cases, it is advisable to 
remove these components from the 
data before the PSD is estimated. 
Without taking this precaution, the 
biasing and strong sidelobes produced 
could easily obscure weaker compo- 
nents. Whenever physically expected, 
the DC component of a signal can 
usually be removed by subtracting the 
sampled data mean, 

1 N ~' 

from each data sample to produce a 
"purely AC" data sequence 

Xq — X, X i — X, . . . , X^ _ i — X • 

ZERO-PADDING THE FFT 

An interesting property of the FFT 
is that simply adding zeros after a 
windowcd-data-samples sequence x,, 
X., x N , to create a longer record x tv 



Listing 1— This program estimates the power spectral distribution ot a complex data sequence using three 
dilterent approaches: the zero-padded FFT, the Welch periodogram method, and Marple's autoregressive 
method. The subroutines listed were adapted from those presented in SL Marple's Digital Spectral Analysis 
with Applications, Englewood cutis, NJ: Prentice Hall, 10S7. 77k. subroutines wm translated to njn under 
QuickBasic 4.5. 



DEFDBL A-Z 



'Use double precision 



COMPLEX ARRAYS 
REM SDYNAMIC: 

DIM ar(l),ai(l).dr(l),di(l).rr(l).ri(l).cr(l),ci(l).xr(l),xi(l) 
DIM w(l).xfftr(l).xffti<l),psd<l).win(l),xsegr(l).xsegi(l) 

CLS 

INPUT "Please input name of data file : ". filename! 
INPUT "Is data complex [y/n] ? ". complexS 
INPUT "Sampling period [seconds] ? ", L 

INPUT "Periodogram number of samples per segment ? ". nsampl 
INPUT "Periodogram number of samples shift ? ". nshiftl 
INPIIl "Auto regressive mnriol order ? ". ip 

' Determine length of data record 
n = 

OPEN filename* FOR INPUT AS #1 
WHILE NOT E0F(1) 
n - n + 1 

IF (complex! - "j" OR complexS - "Y" ) THEN 

INPUT #1. xr(l). xi(l) 
ELSE 

INPUT #1. xr(l) 
END IF 
WEND 
CLOSE #1 



f 



(continued) 
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Listing 1 — continued 



REDIM xr(n). xi(n) ' Redimension data array 

' Read data into array 

OPEN filenames FOR INPUT AS #1 

FOR k - 1 TO n 

IF (complex! - "y" OR complex! = "Y" ) THEN 
INPUT #1. xr(k). xKk) 

ELSE 

INPUT #1, xr(k) 
END IF 
NEXT k 
CLOSE #1 



splay screen on EGA mode 640x350 with 16 colors 




300) - (562, 300) 
250)-(562. 250) 
200) (562. 20(1). 3 
150} - ( 562 . 150), 3 
100 ) - ( 562 . 100). 
50)-(562. 50). 3. 
300)-(45. 50). 3. 
300)1562. 50). 3 



3. 



&H8888 
&H8888 
HIIR8RR 
&H8888 
SH8888 

8.H8888 

&H8888 

, SIIR888 



' Draw di 
SCREEN 9. 
CLS 

LINE (45. 
LINE (45. 
LINE (45. 
LINE (45. 
LINE (45. 
LINE (45. 
LINE (45, 
LINE (562. 

LINE ((562-45J/2+45, 300 ) - ( ( 562 - 45 ) /2+45 . 50), 5. 

LOCATE 22. 2: PRINT "-100"; 

LOCATE 18. 2: PRINT "-80": 

LOCATE 15. 2: PRINT "-60": 

LOCATE 11. 2: PRINT "-40": 

LOCATE 8, 2: PRINT "-20"; 

LOCATE 4, 2: PRINT "0 dB" ; 

LOCATE 3, 28: PRINT "RELATIVE POWER SPECTRUM": 
IF (complex! = "y" OR complex! - "Y" ) THEN 
npsd - 512 

LOCATE 23. 6: PRINT "0.5"; 
LOCAIE 23. 39: PR1NI "0": 
ELSE 

npsd - 1024 

LOCATE 23. 6: PRINT "0"; 
LOCATE 23. 37: PRINT "0.25"; 
END IF 

LOCATE 23.71: PR1NT"0.5"; 

LOCATE 24.18: PRINT"FRACTION OF SAMPLING FREQ. 
LOCATE 1.4: PRINT"PSD ESTIMATORS : "; 
COLOR 11: PRINT"Zero-Padded FFT "; 
COLOR 12: PRINT" Welch's Method ": 
COLOR 10: PRINT"Marple's Method ": 



&H8888 



Fs- 



1/t: 



[Hz]" 



' Compute zero-padded FFT 
nshift - 1 
nsamp = n 
wintype = 1 
GOSUB periodogram 
col - 11: GOSUB plot 



'Set the periodogram for a single segment. & 
'Use rectangular window, in order to 
'Compute zero-padded FFT through periodogram 
'Plot results in light blue 



' Estimate the PSD through Welch's averaged periodogram method 

nshift - nshiftl 'Periodogram num of sample shift between segs 

nsamp - nsampl 'Periodogram num of samples per segment 

wintype - 'Apply Hamming window 

GOSUB periodogram 'Estimate PSD through Welch's method 

col - 12: GOS1IR plot. 'Pint results in light rod 

' Estimate the PSD through Marple's method 

GOSUB marplepsd 'Estimate PSD through Marple's AR method 

col - 10: GOSUB plot 'Plot results in light green 
GOTO progend 

a r p 1 epsd : 

Subroutine to estimate the power spectral distribution of a 
data sequence by Marple's method. This subroutine first solves 
Marple's equations for the estimation of complex autoregressi ve 
coefficients from complex data. Then, it evaluates the transfer 

(continued) 
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x„ x N _,, 0, before performi 
the FFT onuses the FFT to Interpol 
transform values between the N 
original transform values. This 
process, called zero padding, is oft< 
mistakenly thought of as a trick to 
improve the inherent resolution of 
FFT. Zero padding, however, provic 
a much smoother PSD and helps 
annul ambiguities regarding the 
power and location of peaks that m 
be scalloped by the nonzero-, 
FFT. 



CLASSICAL METHODS 

As mentioned before, a commo( 
mistake is to assume that the solutii 
to Equation |2|, the so-called periodc 
gram, is a reliable estimate of PSD. 
Actual proof of this is beyond the 
scope of this article. But, it has been 
demonstrated that regardless of how 
large N is (the number of available 
data samples], the statistical varianc 
of the estimated periodogram spec- 
trum does not .tend to zero. This 
statistical inconsistency is responsib 
for the lack of reliability of the 
periodogram as a spectral estimator. 

The solution to this problem is 
simple, however. If a number of 
periodograms are computed for 
different segments of a data record, 
their average results in a PSD estima 
with good statistical consistency. 
Based on this, Welch [1] proposed a 
simple method to determine the 
average of a number of periodograms 
computed by overlapping segments o 
the available data record. 

Welch's PSD estimate P\f) of M 
data samples is the average of K 
periodograms P(f) of N points each: 



where P(f)'s are obtained by applying 
Equation (2) on appropriately weight! 
data. 

It is obvious that, if the original 
M-point data record is divided into 
segments of N points each, with a sh 
of S samples between adjacent seg- 
ments, the number of periodograms 
that can be averaged is: 

K-Integer(M^) 
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Listing 1 — continued 



' function of the estimated AR system by using the FFT. 

Input Parameters: 
n - Number of data samples (integer) 
ip - Order of linear prediction model (integer) 
xr.xi ■ Array of complex data samps xr(l).xid) - xr(n).xKn) 
npsd - Power spectral distribution length 

Intermediate Parameters: 
p - Real linear prediction variance at order ip 
ar.ai • Array of complex linear prediction coefficients 

' Output Parameters: 

psd - Array containing real power spectral distribution, 
with a maximum power of psdmax 

REDIM ar(ip). ai(ip). dr(ip + 1), ditip + 1). rr(ip). ri(ip) 
REDIM crtip + 1). ci(ip + 1) 

rl - 

FOR k - 2 TO n • 1 

rl - rl + 2 * <xr(k) » 2 + xi(k) « 2) 
NEXT k 

r2 - xr(l) » 2 + xi(l) " 2 
r3 - xr(n) • 2 + xi (n) " 2 
r4 - 1 / (rl + 2 * (r2 + r3)> 
p - rl + r2 + r3 
delta - 1 - r2 * r4 
qamma - 1 - r3 * r4 

Umbdar - r4 * (xr(l) * xr(n) xi(l) * xi(n)) 
lambdai - -rA * (xr(l) * xi(n) - xr(n) * xi ( 1 ) ) 
cr(l) - xr(n) * r4 : ci(l) = xi(n) * r4 
dr(l) - r4 * xr(l): did) - -r4 * xi(l) 
m - 

II ip li 1 1 1 1 M 
p = ( .5 * rl + r2 + r3) / n 

LOCATE 1. 1: BEEP: PRINT "ERROR: Zero AR model order ! 
GOTO progend 
END IF 

Main loop of Marple's Modified Covariance algorithm 

ma rpl el oop : 

n - m + 1 

savelr - 

saveli - 

FOR k = m + 1 TO n 

savelr - savelr + xr(k) * xr(k - m) + xi(k) * xi(k - m) 
saveli - saveli - xr(k) * xi(k - m) + xi(k) * xr(k - m) 

NEXT k 

savelr - 2 * savelr: saveli - 2 * saveli 
rr(m) - savelr: ri(m) = saveli 

thetar - xr( n )*dr( 1 ) -xi < n )*di ( 1 ) : thetai - xr< n )*di ( 1 )+dr( 1 )*xi (n) 
psir = xr(n)*cr(l)-xi(n)*ci(l>: psii - xr(n)*ci ( 1 )+cr( 1 )*xi ( n) 
xir - xr(l)*dr(l)+xi(l)*di(l): xii - xr( 1 )*di ( 1 ) -xi ( 1 )*dr(l ) 

II m <> 1 IIILN 
FOR k - 1 TO m - 1 

thetar - thetar + xr(n-k) * dr(k+l) - xi(n-k) * di(k+l) 
thetai - thetai + xr(n-k) * di(k+l) + xi(n-k) * dr(k+l) 
psir - psir + xr(n-k) * cr(k+l) - xi(n-k) * ci(k+l) 
|)', i i \y. i i i xi (n k ) « i. i ( k i 1 ) i xi (n k ) ' cr(**1 ) 
xir - xir i xr(kH) * dr(kU) * xi(k+l) * di(k+l) 
xii - xii + xr(k+l) * di(k+l) - xi(k+l) * dr(k+l) 
rr( k) - rr( k ) (xr(n+l in) A xr( n+1 m+k )+x i (nH in) *x i ( n+1 -ra+k) ) 
ri(k) - ri (k)-(-xr(n+ 1 -m)*xi (n+1 -m+k)+xi ( n+1 -m)*xr(n+l-m+k) ) 
rr(k) - rr( k) - (xr(m)*xr(m-k)+xi (m)*xi (m-k) ) 
ri(k) = ri ( k) - (xr(m)*xi (m-k) -xi (m)*xr(m-k) ) 
savelr - savelr+rr( k)*ar (m-k )+ri ( k )*ai (m- k) 
saveli = saveli+rr(k)*ai (m-k) ri (k)*ar(m-k) 
NEXT k 

END IF 

(continued) 



HIGH-RESOLUTION METHODS 

The main limitation of FFT-based 
methods is restricted spectral resolu- 
tion. The highest inherent spectral 
resolution (in Hz) possible with the 
FFT is approximately equal to the 
reciprocal of the time interval (in 
seconds) over which data for the FFT is 
acquired. This limitation, which is 
further complicated by leakage and 
the picket-fence effect, is most 
noticeable when analyzing short data 
records. 

It is important to note that short 
data records not only result because of 
the lack of data (such as when sam- 
pling a short transient at a rate barely 
enough to satisfy Nyquist's criterion), 
but also from data sampled from a 
process which slowly varies with time. 

For example, by analyzing the 
vibrations picked up from an oil-well 
drill, the operator can monitor the 
buildup of resonance in the long pipe 
that carries the torque to the drill bit, 
avoiding costly damages to the 
equipment [2]. Although a continuous 
signal from the vibration transducers 
is available for sampling, the vibra- 
tions on the drill assembly change 
rapidly, resulting in a limited number 
of data samples which represent each 
state of the drill bit. It is here where 
high-resolution estimates would be 
desirable, even though the data 
available is limited. 

A number of so-called high- 
resolution spectral estimators have 
been proposed. These alternative 
methods do not assume, as the FFT 
does, that the signal outside of the 
observation window is merely a 
periodic replica of that observed. 
Instead, for instance, the parametric 
estimator relies on the selection of a 
model, which suitably represents the i 
process generating the signal, to 
capture the true characteristics of the 
data outside of the window. By 
determining the model's parameters, 
the theoretical PSD, implied by the 
model, can be calculated and should 
represent the signal's PSD. 

Many signals encountered in real- 
world applications are well approxi- 
mated by a rational transfer function 
model. For example, human speech 
can be characterized by the resonances 
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of the vocal tract that generates it. In 
turn, these resonances are well 
represented hy the poles of a digital 
filter. Parameters for the filter can he 
estimated so that the filter could turn 
a white-noise input into a signal of 
interest. From the filter's transfer 
function, we could easily estimate the 
PSD of the signal. 

Various kinds of filter structures 
exist and are often classified according 
to the type of transfer function they 
implement. An all-pole filter is called 
an autoregressive |AR) model, an all- 
zero filter is a moving-average (MA) 
model, and the general case of a pole- 
zero filter is called an auLoregressive- 
moving-average (ARMA) model. With 
the last example, the model best suited 
for speech is then an AR model. 

Although high-resolution estima- 
tors have heen implemented for all 
these models, AR-based estimators are 
the most popular because many 
computationally efficient algorithms 
are available. A well-behaved set of 
equations to determine the AR 
parameters with a computationally 
efficient algorithm has been intro- 
duced by Marple [3]. 

In the model of Figure 2, the AR- 
filter coefficients a,, a v a are 
estimated by Marple's algorithm based 
on the input data samples x. t 
x N The model assumes that a white- 
noise source drives the filter in which 
the output is regressed through a chain 
of delay elements zr' from which p 
taps feed the AR coefficients. The 
system's transfer function can then be 
computed efficiently through the FFT, 
resulting in an estimate of the signal's 
PSD. 

The performance of Marple's 
estimator is startling. Figure 3a 
presents three spectral estimates 
obtained from a short 64-point com- 
plex-tcst-data set suggested by Marple. 
Estimates obtained through the zero- 
padded FFT periodogram, Welch's 
averaged periodogram, and Marple's 
method can be compared to the 
theoretical spectrum of Figure 3b. 
Only positive-frequency PSD esti- 
mates are shown for clarity. 

Notice that the closely spaced 
components cannot be resolved by 
either of the classical methods, but 



Listing 1— cmtintiad 



-savoli / p 
2) 



tit " • savel r / p: cl 
ar(m) = clr: ai (m) - cli 
p - p * (1 - clr A 2 - cli 
IF m <> 1 THEN 

FOR k - 1 TO m / 2 
mk - m - k 

savelr - ar(k): saveli - ai(k) 
auxr = savelr + clr * ar(mk) + cli * 
ai(k) - saveli - clr * ai(mk) + cli 
ar(k) - auxr 
IF k <> mk THEN 

ar(mk) - ar(mk) 

ai(mk) - ai(mk) 
END IF 



ai(mk) 
' ar(tnk) 



+ clr 
clr 



* savelr + cli * 

* saveli + cli * 



saveli 
savelr 



NEXT k 
END IF 

IF m - ip THEN 

p - .5 * p / (n - m) 

GOTO arpsd 
END IF 

Time update of CD vectors and GAMMA , DELTA , LAMBDA scalars 
rl = 1 / (delta * gamma - lambdar A 2 - lambdai • 2) 
clr - (thetar * lambdar + thetai * lambdai + psir * delta) * rl 
cli = (-thetar * lambdai + thetai * lambdar + psii * delta) * rl 
c2r - (psir * lambdar - psii * lambdai + thetar * gamma) * rl 
c21 - (psir * lambdai + psii * lambdar + thetai * gamma) * rl 
c3r - (xir * lambdar + xii * lambdai + thetar * delta) * rl 



c3i 



(-xir 



c4r - (thetar 



lambdai + xii * lambdar + thetai * 
1 ambdar 



thetai 



c4i - (thetar * lambdai + thetai 
FOR k - 1 TO (m - 1) / 2 + 1 

mk - m + 1 - k 

savelr - cr( k) : 

save2r - dr(k): 

save3r - cr(mk) : 

save4r - dr(mk) : 



1 ambdai 
1 ambdar 



+ xi r 
-t xii 



delta) * 

* gamma) 

* gamma) 



rl 



rl 
rl 



saveli - -ci(k) 
save2i - -di(k) 
save3i - -ci(mk) 
save4i - -di(mk) 
cr(k) - cr(k)+clr*save3r-cli*save3i+c2r*save4r-c2i*save41 
ci(k) - ci (k)+clr*save3i+cli*save3r+c2r*save4i+c21*save4r 
dr(k) - dr(k)+c3r*save3r-c3i*save3i+c4r*save4r-c41*save4i 
di(k) - di (k)+c3r*save3i+c3i*save3r+c4r*save4i+c4i*save4r 
IF k <> mk THEN 

cr(mk) - cr(mk)+clr*savelr-cli*saveli+c2r*save2r-c2i*save2i 
ci(mk) - ci (mk)+cli*savelr+clr*saveli+c2r*save2i+c2i*save2r 
dr(mk) - dr(mk)+c3r*savelr-c3i*saveli+c4r*save2r-c4i*save2i 
di(mk) - di (mk)+c3r*saveli+c3i*savelr+c4r*save2i+c4i*save2r 

END IF 
NEXT k 

r2 - psi r " 2 + psi i A 2 
r3 - thetar A 2 + thetai » 2 
r4 = xi r A 2 + xi i "2 
auxr = psir * lambdar - psii * lambdai 
auxi - psir * lambdai + psii * lambdar 
aux2r - auxr * thetar + auxi * thetai 

r5 - gamma - ( r2 * delta + r3 * gamma + 2 * aux2r) * rl 
auxr = thetar * lambdar - thetai * lambdai 
auxi - thetar * lambdai + thetai * lambdar 
aux2r - auxr * xir + auxi * xii 

r2 - delta - ( r3 * delta + r4 * gamma + 2 * aux2r) * rl 

gamma = r5 
delta - r2 

lambdar = 1 ambdar+c3r*psi r+c3i*psi i+c4r*thetar+c4i*thetai 
lambdai = 1 ambdai -c3r*psi 1+c3i *psi r-c4r*thetai+c4i *thetar 
IF p <- THEN GOTO arpsderr 
IF (delta<-0 OR delta>l OR gamma<-0 OR gammaM) THEN GOTO arpsderr 
rl - 1 / p 

rZ - I I (delta * gamma - lambdar A 2 - lambdai A 2) 

efr - xr(m + 1): efi - xi(m + 1) 

ebr - xr(n - m): ebi - xi(n - m) 

FOR k - 1 TO m 

efr - efr + ar(k) * xr(m + 1 - k) - ai(k) * xi(m + 1 
efi - efi + ar(k) * xi(m + 1 - k) + ai(k) * xr(m + 1 



k) 
k) 

(conlinu&nl 
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Listing 1— continued 

ebr - ebr + ar(k) * xr(n - in + k) + a i ( k ) * xt(n ■ in + k) 
ebi - ebi + ar(k) * xi(n ■ m + k) - ai(k) * xr(n - m + k) 
NEXT k 

clr - ebr * pis cli - ebi * rl 
c2r - efr * rl: c2i - -efi * rl 

c3r - (ebr * delta + efr * lambdar - efi * lambdai) * r2 
c3i - (-ebi * delta + efr * lambdai + efi * lambdar) * r2 
auxr - ebr * lambdar • ebi * lambdai 
auxi - ebr * lambdai + ebi * lambdar 

c4r - (efr * gamma + auxr) * r2: c4i - (efi * gamma - auxi) * r2 
FOR k - m TO 1 STEP -1 

savelr - ar(k): saveli - ai(k) 

ar(k) - savelr+c3r*cr(k)-c3i*ci(k)+c4r*dr(k)-c4i*di(k) 
aKk) - saveli+c3r*ci(k)+c3i*cr(k)+c4r*di(k)+c4i*dr(k) 
cr(k + 1) - cr(k) + clr * savelr - cli * saveli 
ci(k + 1) - ci(k) + clr * saveli + cli * savelr 
dr(k + 1) - dr(k) + c2r * savelr - c2i * saveli 
di(k + 1) - di(k) + c2r * saveli + c2i * savelr 
NEXT k 

cr(l) - clr: ci ( 1 ) - cli 
dr(l) - c2r: did) - c2i 
r3 - ebr * 2 + ebi A 2 
r4 - efr " 2 + efi " 2 

auxr - efr * ebr - efi * ebi: auxi - efr * ebi + efi * ebr 

aux2r - auxr * lambdar - auxi * lambdai 

p - p - (r3 * delta + r4 * gamma + 2 * aux2r) * r2 

delta - delta - r4 * rl 

gamma - gamma - r3 * rl 

auxr - efr * ebr - efi * ebi : auxi - efr * ebi + efi * ebr 
lambdar - lambdar + auxr * rl: lambdai - lambdai - auxi * rl 
IF (p>0 AND delta>0 AND delta<=l AND gamma>0 AND gamma<-l) 
THEN GOTO marpleloop 

arpsderr : 

LOCATE 1, 6: BEEP 

PRINT" ERROR: Numerical i 1 1 -condi t ioni ng detected for model order>"; 

m - 1: 
GOTO progend 

arpsd: 'Evaluate the AR model's transfer function 
nfft - npsd 

RED1M xfftr(nfft). xffti(nirt). psd(npsd) 
xfftr(l) - 1: xffti(l) = 
FOR k - 1 TO ip 

xfftrtk + 1) - ar(k): xffti(k + 1) - ai(k) 
NEXT k 

FOR k - ip + 2 TO npsd 

xfftr(k) - 0: xffti(k) - 'Zero-pad up to npsd 
NEXT k 
GOSUB fft 
psdmax = 
FOR k = 1 TO npsd 

psd(k) - p * t / (xfftr(k) • 2 + xffti(k) « 2) 

IF psd(k) > psdmax THEN psdmax - psd(k) 
NEXT k 
Kl IIIKN 

fft: 

' Subroutine to compute the complex FFT of a complex data series. 

' Input Parameter*.: 

nfft- - FFT size 

t • Sample interval in seconds 

xfftr.xllti Array of tiff I complex data samples 

xfftrd).xffti(l) to xfftr(nfft).xfftKnfft) 

Output Parameters: 

xfftr.xffti - nfft complex transform values replace original 

data samples indexed from k-1 to k-nfft. 

representing the frequencies (k-l)/nfft*t [Hz] 
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they appear clearly separated in the 
estimate produced hy Marple's 
method. You may also notice that 
Marple's estimate is "peaky" even for 
the smooth continuous spectral 
components at the far right of the PSD. 

The reason for this peakiness is 
that a purely autoregressive filter 
generates a spectrum based on pure 
resonances. Only through the use of a 
moving-average could these reso- 
nances be damped to produce a 
perfectly smooth spectrum in regions 
where this is necessary. Although this 
limitation of AR-based estimators 
would lead to errors in the actual 
amplitudes of the PSD components, it 
is very well suited for the high- 
resolution detection of periodicities in 
the signal. 

A price must be paid for the 
increase in resolution and, just as you 
may suspect, the computational 
burden of these high-resolution 
methods far exceed that of a simple 
FFT. In addition, like the selection of 
an appropriate window for the classical 
estimators, the rules for selecting an 
appropriate model, parameter estima- 
tion method, and model order are all 
but cast in stone. 

IMPLEMENTING SPECTRAL 
ANALYSIS ALGORITHMS 

Program spectrum. bas presented 
in Listing 1 demonstrates the imple- 
mentation of the spectral estimation 
methods discussed. The program was 
written in QuickBASIC 4.5, but should 
run with little trouble under any other 
BASIC compiler on an IBM PC- 
compatible with EGA/VGA graphics. 
However, BASIC does not support 
complex-number arithmetic, so 
explicit operations have been used in 
which variable names with the suffix i 
represent the real portion of that 
variable, while those with the suffix i 
represent the imaginary portion of the 
same. 

After being defined by the user, 
the program reads a file containing the , 
N-data-point sequence to be analyzed. . 
The data can be either a single column 
of (plain ASCII) samples or two 
columns, one containing the real and 
the other, the imaginary parts of 
complex data samples. The program 
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PC-Based Instruments 
200 MSa/s DIGITAL OSCILLOSCOPE 



HUGE BUFFER 
FAST SAMPLING 
SCOPE AND LOGIC ANALYZER 
C LIBRARY W/SOURCE AVAILABLE 
POWERFUL FRONT PANEL SOFTWARE 



DSO Channels 

2 Ch. up to 100 MSa/s 

or 

1 Ch. at 200 MSa/s 
4K or 64K Samples/Ch 
Cross Trigger with LA 
125 MHz Bandwidth 
Logic Analyzer Channels 
8 Ch. up to 100 MHz 
4K or 64K Samples/Ch 
Cross Trigger with DSO 




$1799 - DSO-28204 (4K) 
$2285 - DSO-28264 (64K) 



Universal Device Programmer 



PAL 

GAL 

EPROM 

EEPROM 

FLASH 

MICRO 

PIC 

etc.. 




$475 



Free software updates on BBS 
Powerful menu driven software 



400 MHz Logic Analyzer 



up to 128 Channels 
up to 400 MHz 
up to 16K Samples/Channel 
Variable Threshold Levels 
8 External Clocks 
■ 16 Level Triggering 
Pattern Generator Option 



$799 . LA12100 (100 MHz, 24 Ch) 
$1299 - LA32200 (200 MHz, 32 Ch) 
$1899 - LA32400 (400 MHz, 32 Ch) 
$2750 - LA64400 (400 MHz, 64 Ch) 




A 



Call (201) 808-8990 

Link Instruments 

%\f 369 Passaic Ave, Suite 100, Fairfield, NJ 07004 fax: 808-8786 
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will estimate the spectrum of the 
input dfrt.1 using throe methods: 

1 1 A single periodogram of the data 
record is obtained by zero padding 
the data up to npsd data points 
[npsd = 512 for complex and 1024 
for real input data) from which the 
squared modulus of the FFT is 
computed. A rectangular window i 
assumed. 

2| Welch's method with a Hamming 
window is applied using the 
number of samples per periodogran 
and the shift specified by the user. 

3) Marple's method is used to estimati 
the PSD of the data using an AR 
model with model order given by 
the user. 

Prior to its display in the output 
screen, PSD is normalized relative to 
its maximum, and transformed to 
decibels. For complex input data, both 
the positive and negative frequency 
sides of the spectrum are plotted. 
Otherwise, only the positive frequenq 
spectrum is presented. 

Because of screen resolution 
limitations, the number of computed 
PSD points for display has been 
limited to 512. If a larger PSD record is 
required, however, npsd can be 
increased to any desired power of 2, 
and a file can be opened to receive the 
PSD-estimate results. 

A few simple demonstrations can 
be set up to compare the performance 
of the methods. First, you may 
generate a data file for a signal consist 
ing of a single sinusoid at Vtf t with 
white noise added to it using the 
program in Listing 2. 

You may vary the signal-to-noise 
ratio by changing the value of the 
noise component's coefficient. As 
well, the frequency of the sinusoidal 
component may be changed by alterinj 
the denominator of the sine argument 
Of course, from Nyquist's theorem, a 
denominator smaller than two pro- 
duces an aliased signal (you may want 
to experiment with the effect that this 
has on the PSD estimate|. 

In addition, the resolving power of 
the estimators can be compared by 
using a signal containing two closely 
separated sinusoidal components. This 
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xffti(kk) 
d l-xffti(kk) 



Listing 1— continued 

RED I M wr(nfft) AS DOUBLE . wi(nfft) AS DOUBLE 

' Set up complex exponential table Tor rn 

nexp - 1 

nt - 2 A nexp 

WHILE nt < nfft 

nexp - nexp + 1 

nt - 2 " nexp 
WEND 

IF nt <> nfft THEN 

LOCATE 1. 4: BEEP : PRINT "Error! nfft is not a power of 2 

GOTO progend 
END IF 

s * 8 * ATN(l) / (nt) 

clr - COS(s): cli - -SIN(s) 

c2r - 1: c2i - 

' Compute complex exponential array 
FOR k - 1 TO nt 

wr(k) - c2r: wi(k) - c2i 

auxr - c2r * clr - cZ i * cli: c2i - c2r * cli + c2 i » clr 
c2r - auxr 
NEXT k 



' Main FFT routine 

mm - 1 

11 - nfft 

FOR k - 1 TO nexp 

nn - 11 / 2 

jj - mm + 1 

FOR i - 1 TO nfft STEP 11 
kk - i + nn 

clr - xfftr(i) + xfftr(kk): cli - xffti(i) + 
xfftr(kk)-xfftr(i J-xfftr(kk): xf f tl ( kk )-xf f ti 
xfftr(i) = clr: xffti(i) - cli 
NEXT i 

IT nn <> 1 THTN 
FOR j - 2 TO nn 

c2r - wr( jj ): c2i - wi ( jj ) 
FOR i = j TO nfft STEP 11 
kk - i + nn 

clr - xfftrd )+xfftr(kk) 
auxr - xf ftr( i ) -xf ftr(kk 
xfftr(kk) - auxr * c2r 
xffti(kk) - auxr * c2i + auxi * c2r 
xfftrd) - clr: xffti(i) - cli 
NEXT i 

jj " jj + mm 
NEXT j 
11 - nn 
mm - mm * 2 
END IF 
NEXT k 

nv2 - nfft / 2 
nml - nfft - 1 
j " 1 

FOR i - 1 10 mill 
II i ■ j Mil N 

Clr - xfftr(j): cli - xffti(j) 
xfftr(j) - xfftrd): xffti(j) - xfftid) 
xfftrd) = clr: xfftid ) - cli 
LNI) II 
k - nv2 
Willi I k < j 
j - j • k 
k - k / 2 
WEND 

3 - j + k 
NEXT 1 

FOR i - 1 TO nfft 

xfftr(i) - xfftr(i) - t: xfftid) - xfftid) * t 
NEXT i 
RETURN 



cli = xfftid 
: auxi = xffti 
auxi * c2i 



)+xffti(kk) 
d)-xffti(kk) 
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NEW Data 

Acquisition 

Catalog 

Covers expanded 
low cost line. 



FREE! 



1994 120 page catalog for PC, VME, 
and Qbus data acquisition. Plus infor- 
mative application notes regarding 
anti-alias filtering, signal condition- 
ing, and more. 



LabVIEW", LabWindows", 
Snap-Master ™, and more 

NEW Low Cost I/O Boards 

NEW Industrial PCs 

NEW Isolated Analog and 
Digital Industrial I/O 



New from the inventors of 
plug-in data acquisition. 

Call, fax, or mail for your 
free copy today. 



ADAC 



American Data Acquisition Corporation 
70 Tower Office Park, Woburn, MA 01801 
Phone: (800) 648-6589 Fax: (617) 938-6553 
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Replace Four 
Conventional PC/104 
Modules with 
One SuperXT'" 
CMF8680 cpu Module " 

Embedded PC/XT Controller with 
Intelligent Power Management 




■ PC/XT compatibility with 286 emulation 

■ 14 MHz, 16-bit 8086 CPU 

■ +5V only; 1 .6W at 1 4.3 MHz, 1W at 7.2 MHz 

■ Intelligent sleep modes, 0.1W in Suspend 

■ ROM-DOS and RTD enhanced BIOS 

■ Compatible with MS-DOS & real-time 
operating systems 

■ 1M bootable Solid State Disk & free software 

■ 4K-bit configuration EEPROM (2K for user) 

■ 2M on-board DRAM 

■ IDE & floppy interlaces 

■ CGA CRT/LCD controller 

■ Two RS-232 ports, one RS-485 port 

■ Parallel, XT keyboard & speaker ports 

■ Oplional X-Y keypad scanning/PCMCIA 
interface 

■ Watchdog timer & real-time clock 

Expand This Or Any PC/104 System 
with the 
CM106 SuperVGA 
Controller utilityModule™ 

■ Mono/color STN & TFT flat panel support 

■ Simultaneous CRT & LCD operation 

■ Resolution to 1024 x 768 pixels 

■ Displays up to 256 colors 



$223 



Speed Product Development with the 
DS8680 Development System 

Your DS8680 includes the CMF8680, CM102 
keypad scanning/PCMCIA, CM104 with 1.8" 
85MB hard drive, CM106 SVGA controller & 
DM5406 12-bit. 100 kHz dataModule™ in an 
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Listing ^-continued 



pprintlogram: 

Subroutine t.o compute averaged periodogram over nseg segment; 
Input Parameters: 
n - Number of data samples 

nshift - Number of samples shift between segments 
nsamp - Number of samples per segment (must be even) 
t - Sample interval in seconds 

xr.xi - Array of complex samples xr(l),xi(l) to xr(n).xi 
wintype - Window type - - Hamming, other - rectangular 



Output Parameters: 

nseg - Number of segments averaged 

psd - Array containing real power spectral distribution 
with a maximum power of psdmax 

REDIM win(nsamp) , xsegr(npsd). xsegi(npsd), psd(npsd). 

xfftr(npsd), xffti(npsd) 
pi 2 - 8 * ATN(l) 

' Compute window 
FOR k = 1 TO nsamp 
IF wintype - THEN 
' Hamming window 

win(k) = 0.538 + 0.462 * C0S( pi 2*( -0. 5+( k- 1 ) /(nsamp- 1 )) ) 
ELSE 

win(k) - 1 'Rectangular window 
END IF 
NEXT k 

' Compute Welch's averaged periodogram applying window 
nseg - (n - nsamp) / nshift + 1 
FOR k - 1 TO nseg 
FOR j - 1 TO nsamp 

index - j + (k - 1) * nshift 
xsegr(j) - xr(index) * win(j) 
xseyi(j) - xi( index) * win(j) 
NEXT j 

FOR j - nsamp + 1 TO npsd 

xsegr( j ) - 0: xsegi ( j ) = 
NEXT j 
nfft - npsd 
FOR j - 1 TO nfft 

xfftr(j) - xsegr(j): xffti(j) - xsegi(j) 
NEXT j 
G0SUB fft 
FOR j - 1 TO npsd 

IF k - 1 THEN 

psd(j) - xfftr(j) » 2 + xffti(j) * 2 

ELSE 

psd(j) - psd(j) + xfftr(j) » 2 + xffti(j) - 2 
END IF 

NEXT j 
NEXT k 
psdmax - 
FOR k = 1 TO npsd 

psd(k) = psrt(k) / (nseg * nsamp) 

IF psd(k) > psdmax THEN psdmax - psd(k) 
NEXT k 
RETURN 



f 



'Zero-pad up to npsd 



plot: 



assuming npsd - 512 



Plot results on graph using color col 
(complex) or npsd - 1024 (real) 
FOR k - 1 TO npsd 

psd(k) = 20*L0G(psd(k)/psdmax)/L0G(10) 
IF psd(k) < -100 THEN psd(k) - -100 
NEXT k 

IF (complex$ - *y" R complex! - "Y" ) THEN 
' Plot PSD for positive frequencies 
FOR k - 2 TO 256 

LI NE(44+k+256, 50-2. 5*psd( k- 1 ) ) - (45+k+256 . 50-2. 5*psd(k)). col 



Normal ize S. xform data 
Clip at -100 dB t 
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Listing \— continued 



NEXT k 

' Plot PSD for neqative frequencies 
FOR k - 2'j(> 10 !>1Z 

LINE(44+k- 256, 50-2. 5*psd( k- 1 ) ) - (45+k-256 , 50-2. 5*psd(k)). col 
NEXT k 
ELSE 

' Plot PSD only for positive frequencies 

FOR k - 2 TO 512 

LINE (44+k. 50-2.5*psd(k-l ))-(45+k, 50- 2 . 5*psd( k ) ) , col 

NEXT k 
END IF 
RETURN 

progend: 

' Wait for any key to be pressed to exit program 
WHILE INKEYS - ": WEND 

END 




Figure 4— Frequency-wavenumber spectra ot ptanewaves oilers (a) spaliotemporal (left) and frequency- 
wavenumber (right) representations ola sinusoidal wave of frequency!, traveling at propagation velocity v,, (b) 
sinusoidal wave ol frequency I, traveling at propagation velocity v,, (c) a sum ol the above. The two-dimensbnal 
PSDs clearly show the component waveform spectra and propagation velocity. 



can be accomplished by adding a 
second component to the program line 
which computus x as, for example: 

x-2(rnd-0..S) + (sin(2n^)) + (»in(2jt- 'J) 

A rule of thumb is to keep the AR 
model order smaller that one-third of 
the number of available data samples 
and to allow for at least twice the 
number of expected spectral compo- 
nents. The code in Listing 1 announces 
an error whenever mathematical ill- 
conditioning is encountered due to 
too-high a model order. However, an 
unreasonably "peaky" spectrum is 
often obtained before ill conditioning 
can be detected. 

AN APPLICATION: 

ARRAY SIGNAL PROCESSING 

The greatest interest in high- 
resolution spectral estimators has 
been generated in the field known as 
array signal processing. Here, a 
number of sensors are placed at 
various locations in space to detect 
traveling waves. 

For example, in seismology, a 
number of sensors capable of detecting 
the shock waves of a tremor or 
earthquake are spread over a certain 
area. As the shock waves travel under 
the sensor array, signals from each 
sensor can be sampled in real time, 
producing a data record which also 
contains information regarding the 
spatial characteristics of the waves 
(because the sensor locations are 
known]. The processing of resulting 
spatiotemporal data is meant to 
estimate the number, vector velocity 
(speed and direction), and wave shape 
of the overlapping traveling waves 
in the presence of interference and 
noise. 

Array signal processing has been 
successfully applied to many fields 
besides seismology (4(. For example, in 
exploration seismogeology, it has been 
used to image limited regions of the 
Earth's interior. In passive sonar, it can 
determine the location and signature 
of submerged sound sources; in radar, 
it acquires targets with very high 
spatial resolution,- and in biomedical 
diagnosis, it tracks weak electrical 
potentials from the brain, nerves, and 
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muscles. Other applications include 
image reconstruction from projections 
such as radio astronomy and medical 
tomography. 

The most common form of 
traveling wave is the planewave. In its 
simplest form, a planewave is a 
sinusoidal wave that not only propa- 
gates through time £, but also through 
space. In the direction of propagation r, 
this wave can be represented by: 



g|t,r| = Asin(27i(ft-X)) 



(6) 



where A is the amplitude of the wave, 
/ is its temporal frequency (Hz = '/.), 
and v is the velocity (in m A or other 
suitable velocity units) at which the 
wave propagates through space. 

If one such simple planewave is 
sampled discretely along time and 
space, we would obtain a record 
similar to that presented in the left 
side of Figure 4a. As you see, at any 
given time the spatial sampling of the 
wave also forms a sinusoid with 
frequency k r The spatial frequency (in 
Vm) of such a simple planewave is 
called the wave number, and is given 
by: 



.1 



[7| 



Its physical meaning indicates that at a 
distance i from the origin, the phase of 
the wave accumulates by 2%kr radians. 
The two-dimensional spectrum of the 
planewave in our example would be an 
impulse S (the spectrum of a sinusoid) 
located in the frequency-wave number 
if-k) plane at ft . Through this kind 
of spectral analysis, we infer the 
components of the waveform and their 
velocity because the slope at which 
the components arc found is equal to 
their propagation velocity or, in this 



".[fi- 



ll 



Adding a second component 
(Figure 4b) with a different frequency 
and propagation velocity to the 
original component, we obtain a 
planewave (Figure 4c| that, regardless 
of its simplicity, can hardly be recog- 
nized in the space-time domain. 



Listing 2— Ms program generates a lile containing data synthesized horn a single sinusoidal signal 
contaminated by random noise. 



pi - 1 . 1 1 1 Wf.? 

OPEN "noise.dat" 

FOR i - 1 TO 256 
x - 2 * (RND - 
PRINT #1, x 

NEXT 1 

CLOSE #1 



TOR OUIPUI AS til 

.5) + (SINC2 * pi * i / 4)) 



However, the two-dimensional 
frequency-wave number spectrum of 
the signal clearly resolves the compo- 
nents and their propagation velocities. 

The two-dimensional spectrum 
can be computed with ease knowing 
that the two-dimensional DFT is 
computed as a sequence of one- 
dimensional DFTs of the columns of 
the data array, followed by a sequence 
of one-dimensional DFTs of the rows 
of this new array, or vice versa. As 
such, the most simple two-dimen- 
sional PSD estimator is implemented 
through the FFT. In practice, however, 
due to the limited number of spatial 
samples (because only a few sensors 



are normally used), the use of high- 
resolution estimators is essential. 

Considering that enough samples 
x,„ x,, x N , can usually be obtained 
from each of the R sensors through 
time, a hybrid two-dimensional 
spectral estimator can be implemented 
by combining the classical and the 
high-resolution spectral estimation 
approaches. As shown in Figure 5, 
using spatiotemporal data g(t, r), 



g(l,r|. 



Xo.R-l X 1,R-1 



X N-1,0 
X N-1,1 



XN-l.R-l 



b) 



g(t.O 

(linear array data) 



g(t.O) FFT > Guff.O) 



9 0.1) - 
g (t.R-1) - 



G„(M) 

G M (f.R-D 



G*(f.r) 
(complex) 



g(U) 



J 



2-D PSD 



G(l.k) 





Figure 5— These images illustrate a hybrid two-dimensional spectral estimator. Spatiotemporal data (a) is 
transformed along time-domain into an intermediate array (b) through the application of a windowed FFT to each 
and every row ot the original data. Applying an AR-PSD estimator to every column ol the intermediate array 
completes the two-dimensional PSD estimation process. 



Issue #52 November 1 994 The Computer Applications Journal 




tins [usee] 



HU Territory 



FILE NUUB.HAP 






Tendon position 




70 Cm] 


Fiber type 




Constant 


Nunber of fibers 




100 


Snatial sanding 




2.34 Inn] 



Conduction velocity 




Figure 6— Frequency-wavenumber spectral estimation has been applied to the analysis ol the potentials recorded 
from a muscle twitch. In (a), the complex spatiotemporal waveform has been analyzed to show information about the 
conduction velocity, origin, and location ol the component potentials; tb) shows a magnified view of the spatiotempo- 
ral data. 



an intermediate transform G lltl [f, r) is 
computed by applying the FFT along 
each row (time domain) of appropri- 
ately weighted data. The two-dimen- 
sional spectral estimate g[f, k) is then 
completed by obtaining the AR-PSU of 
each column of complex numbers in 
the intermediate transform. 

In the more general case, using an 
array of sensors spread out over an area 
with a planewave traveling in any 
direction under the array, a three- 
dimensional hybrid spectral estimator 
determines not only the wave's 
components and its velocities, but also 
each component's bearing. 



For example, tiny electrical 
potentials can be picked up from 
muscle fibers using electrodes attached 
to the skin. These potentials are 
caused by pulses (action potentials! 
that travel down every muscle fiber 
causing the contraction of muscles. 
The conduction velocity as well as the 
origin of these potentials enclose a 
wealth of information which can be 
used as an aid in the early diagnosis of 
nerve and muscle diseases. The large 
number of convoluted signals and the 
very small differences between their 
waveforms makes it impossible to 
determine this information from 



spatiotemporal data (Figure 6b). 
However, a complete analysis (Figure 
6a) is possible through the use of 
multidimensional spectral estimates. 

IN CONCLUSION 

Of course, the BASIC program 
listed here may be too slow to cope 
with most real-time applications, but 
implementing both classical and high 
resolution methods on DSP is a 
relatively easy task. First of all, 
modern DSP chips are specifically 
designed to perform the convolution, 
vector arithmetic, and FFT operations 
in a minimal number of clock cycles. 
In addition, optimized subroutines 
implementing the most popular high- 
resolution algorithms are available 
often in the public domain. 

Multidimensional PSD estimatioi 
has a very high intrinsic parallelism 
because spectral estimates are taken 
independently for every dimension 
and, as such, can be solved efficiently 
in parallel. In other Words, since tasks 
in array-signal processing require 
specific operations to be performed on 
innumerable data blocks, a parallel 
system exploits the full power of a 
number of processors working con- 
comitantly on different portions of thi 
data to solve the larger problem 

High-power computational 
engines (e.g., Intel's i860) and floating- 
point DSPs (e.g., Texas Instruments' 
TMS320C30 and the AT&T DSP32) 
possess the raw floating-point perfor- j 
mance necessary to efficiently imple- 
ment the relevant algorithms. Unfor- 
tunately, however, these chips do not 
present the flexibility required to 
implement multiprocessor architec- 
tures which can optimally exploit 
intrinsic parallelism. Moreover, 
parallel DSP systems using these chip 
would most likely encounter serious 
communication bottlenecks imposed 
by their classical bus-based architec- 
tures. In these cases, RISC chips, such 
as the Transputer family, or DSP 
chips, such as the TM5320C4Q, which 
are designed for parallelism, display 
the full power of a scalable and very 
flexible architecture. 

I have tried to show you that 
spectral analysis is a very convenient 
tool that serves a number of engineer- 



38 



Issue #52 November 1994 The Computer Applications Journal 



ing applications. Moreover, with 
today's PCs, you have the power to 
implement modern PSD estimation 
algorithms with sufficient efficiency 
for experimenting and even for some 
real applications. With the enhanced 
capabilities of DSP chips, PCs with 
DSP coprocessors and laboratory 
spectrum analyzers with embedded 
DSPs become truly powerful and 
useful instruments. 

However, as you understand by 
now, obtaining good spectral estimates 
is not only a matter of blindly applying 
the algorithm anil watching the screen. 
Rather, knowledge about the spectral 
estimation methods and empirical 
experience of their use are of foremost 
importance in obtaining consistent 
results. H 

David Prulchi has a Ph.D. in Biomedi- 
cal Engineering from Tel-Aviv Univer- 
sity. He is an engineering specialist at 
Intermedics, and his main Ra)D 
interest is biomedical signal process- 
ing in implantable devices. He may be 
reached at davidiy@mails.imed.com. 
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Data Genie offers a full line of test & measure- 
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